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ABSTRACT 

The  University  of  Florida,  Electronic  Communications  Laboratory  (ECL),  as  part  of  a  project  with  the  Air  Force  Research 
Laboratory  (AFRL),  is  investigating  the  utility  of  ground  penetrating  radar  (GPR)  for  airfield  evaluations.  It  is  expected 
that  the  GPR  and  the  results  of  current  research  will  assist  the  Air  Force  in  assessing  airfields  with  less  coring  and  allow 
more  rapid  and  accurate  airfield  evaluations  by  providing  continuous  estimates  of  the  pavement  and  subsurface 
conditions.  AFRL  has  conducted  airfield  measurements  using  a  high-resolution  GPR  mounted  on  a  test  van.  Using  time 
domain  data  from  a  pair  of  ground-coupled  antennas,  one  in  bistatic  mode  and  the  other  in  monostatic  mode,  the  thickness 
of  pavement  and  other  subsurface  layers  can  be  estimated.  In  order  for  GPR  to  be  useful  for  airfield  evaluations,  practical 
considerations  such  as  automated  data  processing  capability  and  efficient  data  analysis  procedures  have  also  been 
addressed.  This  paper  will  describe  novel  approaches  and  signal  processing  techniques  developed  for  detecting  layer 
interface  boundaries  and  estimating  the  layer  properties  using  GPR.  In  addition,  the  software  tools  developed  to  facilitate 
airfield  evaluations  will  be  discussed.  Finally,  surface  layer  thickness  estimates  will  be  compared  to  measured  cores  from 
a  recent  airfield  assessment. 

Keywords:  ground  penetrating  radar,  radar  signal  processing,  radar  detection 

1.  BACKGROUND 

Currently,  the  Air  Force  Civil  Engineering  Support  Agency  (AFCESA)  and  other  Air  Force  personnel  perform  airfield 
studies  ranging  from  basic  field  investigations  to  detailed  and  complex  airfield  pavement  evaluations  required  by  the  U.S. 
Air  Force  Airfield  Evaluation  Program  .  The  purpose  of  this  program  is  to  determine  the  condition  of  all  airfields  used  by 
present  or  potential  Air  Force  missions.  Current  pavement  and  subsurface  field  assessment  methods  include  core  drilling, 
visual  observations,  and  obtaining  measurements  using  a  heavy  weight  deflectometer  (HWD)  and  a  soil  penetration 
device  such  as  a  dynamic  cone  penetrometer  (DCP).  These  methods  are  time  consuming  and  can  fail  to  provide  the 
sampling  resolution  needed  to  correctly  characterize  the  surrounding  pavement  and  subsurface  conditions,  track  changes 
in  the  pavement  and  subsurface,  or  identify  localized  problem  areas  underneath  the  pavement.  GPR  has  the  potential  to 
improve  airfield  evaluations  by  allowing  more  continuous  and  rapid  data  collections  than  the  current  evaluation  methods. 
When  combined  with  current  evaluation  techniques  GPR  will  enable  better  decisions  and  help  provide  a  more  accurate 
and  complete  analysis  of  the  pavement  subsurface. 

GPR  is  a  short-range  sensor  that  is  used  to  measure  the  electromagnetic  reflections  from  an  ultra-wideband  radar  impulse. 
The  observed  amplitude  changes  in  the  received  time  response  waveform  -are  due  to  the  abrupt  changes  in  the  electrical 
properties  of  a  layered  medium  or  reflections  from  buried  objects2.  Antennas  used  for  GPR  typically  have  a  beamwidth  of 
90  degrees  or  more.  This  can  make  the  position  of  an  object  or  layer  very  difficult  to  determine  when  considering  only  a 
single  radar  return.  Therefore,  GPR  data  is  typically  collected  along  a  straight  path  in  evenly  spaced  intervals  using  a 
distance  measuring  instrument  (usually  a  survey  wheel)  as  a  trigger. 
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The  ECL  has  worked  with  the  Florida  Department  of  Transportation  (FDOT)  for  several  years  to  improve  the  usability 
and  accuracy  of  GPR  for  roadway  analysis3*4.  The  potential  utility  of  GPR  for  pavement  measurements  has  been  only 
partially  realized  for  airfield  applications.  Several  obstacles  have  prevented  wider  GPR  utilization,  including  the  large 
amounts  of  data  that  need  to  be  collected  by  the  GPR  to  map  an  entire  airfield,  the  significant  amount  of  operator 
interaction  and  interpretation  required,  and  the  vast  amount  of  signal  processing  that  is  required.  These  obstacles  have 
tended  to  make  GPR  under-utilized,  even  in  those  cases  where  adequate  measurement  accuracies  have  been  possible.  The 
goal  of  the  current  research  is  to  improve  the  usability  and  accuracy  of  GPR  for  airfield  assessments. 


2.  GPR  ANALYSIS  SOFTWARE 

As  part  of  the  previous  work  for  the  FDOT,  a  comprehensive  GPR  analysis  software  tool  was  developed  at  the  ECL.  The 
GPR  software  was  designed  and  constructed  to  improve  GPR  evaluation  by  providing  better  interaction  with  the  data  and 
by  increasing  layer  property  estimation  accuracy  to  reduce  subjective  operator  analysis.  The  software  is  driven  by  a 
Graphical  User  Interface  (GUI)  and  provides  a  multitude  of  functions  for  GPR  data  processing,  visualization,  and  data 
management.  For  airfield  assessments,  modifications  were  made  to  the  GPR  analysis  software  to  accommodate  the 
AFRL’s  current  and  potential  GPR  hardware.  Also,  software  tools  were  added  for  batch  processing  and  summary  analysis 
report  creation  for  entire  airfields.  Finally,  import  and  display  capabilities  were  implemented  for  core  log  files  and  data 
from  HWD  and  DCP  sensors.  Figure  1  below  shows  the  main  interface  of  the  GPR  analysis  software. 


Figure  1.  GPR  analysis  software 
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3.  GPR  SIGNAL  PROCESSING  OVERVIEW 


The  main  processing  objective  for  GPR  data  is  to  determine  the  properties  of  subsurface  layers  (asphalt,  concrete,  and  the 
subgrade).  Autonomous  determination  of  these  semi-continuous  subsurface  interfaces  is  a  difficult  task  that  greatly 
impacts  GPR  analysis  efficiency.  Manual  selection  of  interface  returns  is  labor  intensive,  emphasizing  the  need  for 
algorithms  that  allow  the  bulk  of  this  work  to  be  handled  by  a  computer.  GPR  signal  processing  attempts  to  remove  much 
of  the  intensive  human  interaction  and  experience  that  is  required.  This  will  reduce  the  subjectivity  of  operator  analysis 
and  minimize  processing  cost  by  lowering  the  required  operator  training  and  experience.  Several  processing  aspects  are 
used  to  simplify  GPR  analysis.  First,  automated  processing  of  GPR  data  will  attempt  to  detect  and  segment  layer 
interfaces.  Next,  properties  of  the  layers  will  be  computed,  including  layer  thickness,  relative  dielectric  permittivity  (e^, 
and  waveform  propagation  velocity.  Finally,  processing  results  can  be  displayed  and  summarized  using  the  GPR  analysis 
software. 

Techniques  for  layer  property  determination  are  subject  to  the  physics  of  the  GPR  phenomena.  Layer  thickness,  for 
example,  is  estimated  by  measuring  the  delay  and  waveform  velocity  to  each  layer.  Velocity  is,  therefore,  a  critical 
property  which  may  be  determined  in  different  ways.  Air-launched  GPR  systems  can  use  the  ratio  of  the  return  amplitudes 
from  successive  layers  to  estimate  the  reflection  coefficient  (p),  e,.  and  velocity3*5.  This  calculation  may,  at  times,  be 
subject  to  recursive  errors  that  originate  with  surface  amplitude  errors  caused  by  surface  roughness4,6.  An  alternative 
approach  has  been  developed  that  uses  the  time-delay  differences  between  two  ground-coupled  antennas.  This  method  has 
the  potential  for  improved  velocity  and  er  measurements,  which  will  result  in  more  accurate  thickness  calculations. 
However,  processing  data  from  the  dual  propagation  path  system  presents  its  own  unique  challenges. 

GPR  data  processing  is  done  using  a  set  of  sequentially  applied  algorithms  which  automate  the  computation  of  the 
subsurface  properties  needed.  Many  of  the  algorithms  operate  as  general  purpose  functional  processing  blocks  that  work 
independent  of  radar  type  or  configuration.  Other  processing  steps  require  algorithms  that  are  tailored  for  specific  systems 
or  configurations.  The  processing  flow,  as  indicated  in  Figure  2,  has  been  designed  and  implemented  in  a  modular  fashion 
to  facilitate  algorithm  development  and  performance  comparisons  using  the  GPR  analysis  software.  Major  processing 
steps  include  detecting  radar  returns,  prescreening  detects  from  layers,  segmenting  detects  into  interfaces,  improving 
interface  continuity,  and  calculation  of  layer  parameters. 
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Figure  2.  GPR  signal  processing  flow  diagram 

4.  DETECTION  AND  TRACKING 


4.1  Filter  selection 

Matched  filtering  is  a  standard  approach  to  conventional  radar 'signal  processing  and  is  mathematically  identical  to  a 
correlation  receiver7.  The  matched-filter  waveform  is  convolved  with  the  received  GPR  waveform  to  maximize  the  signal- 
to-noise  ratio  (SNR).  Determination  of  the  appropriate  matched  filter  to  use  for  detection  may  be  accomplished  in 
different  ways  depending  on  the  system  configuration.  The  simplest  way  of  acquiring  a  matched  filter  is  to  measure  the 
transmitted  pulse.  The  measured  response  from  an  isolated  interface  is  time-inversed  to  create  a  matched-filter  waveform. 
Measurement  of  an  isolated  interface  response  can  be  easily  accomplished  with  an  air-launched  system.  A  metal  plate 
may  be  placed  on  the  pavement  surface  under  the  elevated  air-launched  antenna.  The  metal  plate  reflects  all  transmitted 
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energy  and  provides  a  measured  return  that  is  free  from  interference.  Ground-coupled  systems  make  field  measurements 
of  the  transmitted  waveform  impractical  due  to  the  inability  to  easily  isolate  an  interface  return  from  the  coupling 
response  and  interfering  returns  within  the  ground.  To  overcome  this  difficulty,  a  GPR  testbed  was  constructed  that  allows 
the  response  from  a  metal  plate  through  a  homogenous  media  to  be  collected.  Isolated  interface  responses  were  obtained 
using  ten  inches  of  concrete  for  each  of  AFRL’s  ground-coupled  antennas.  When  processing  a  particular  data  set,  the 
appropriate  matched-filter  waveform  is  selected  based  on  the  antenna  identification  listed  in  the  GPR  data  header. 

4.2  Detection 

The  received  GPR  signal  resulting  from  a  transmitted  pulse,  x(t),  emerging  from  N-Iayers  of  isotropic  and  homogeneous 
media,  such  as  surface  and  subsurface  layers,  can  be  written  as 

*r(  0  =  2f=  0Atx(t  - 1{)  +  c(t)  +  n(t)  (1) 

where  Af  represents  the  amplitude  of  the  reflected  pulse  from  layer  i  at  time  tif  c(t)  represents  a  noise  component  due  to 
clutter  responses,  and  n(t)  is  a  random  noise  process.  Detection  can  be  viewed  as  the  process  of  finding  the  A,-  and  /,*  from 
xr(t)  in  (1).  Subsurface  layers  are  defined  as  extended  regions  within  the  ground  where  abrupt  changes  in  er  occur.  The 
change  in  ^  results  in  a  large  impulse  response  at  the  interface  boundary.  Layer  interface  detection  attempts  to  obtain  the 
location  and  amplitude  of  the  ground  response  components  associated  with  these  subsurface  layers.  When  GPR  data  is 
processed,  the  number  of  desired  interfaces  to  detect  is  passed  as  a  parameter  to  the  algorithm. 

An  iterative  matched-filter  technique  has  been  used  to  find  the  peak  locations  corresponding  to  interfaces  by  convolution 
of  an  antenna-specific  matched-filter  waveform  with  the  signal  response  xr(t).  At  those  time  sample  points  for  which  a 
peak  in  the  convolution  occurs,  the  sample  point  and  the  sign  of  the  convolution  are  noted.  A  replica  of  the  matched-filter 
waveform  is  scaled  to  the  amplitude  of  the  signal  response,  shifted  in  time  to  the  sample  point  corresponding  to  peak 
correlation  and  subtracted  from  the  signal  response.  This  convolution  and  replica  subtraction  process  is  repeated  until 
peaks  corresponding  to  all  possible  interfaces  have  been  found. 

After  detection,  the  relative  values  of  the  matched-filter  correlation  to  the  change  in  Er  are  disproportionate  due  to  the 
progressive  signal  attenuation  through  the  ground  as  the  distance  from  the  GPR  antenna  increases.  Ideally,  a  detection 
confidence  value  will  be  similar  for  detects  that  correspond  to  similar  ^  changes  regardless  of  distance  from  the  GPR’s 
antenna.  To  overcome  this  problem,  a  model  of  the  media’s  attenuation  characteristics  needs  to  be  developed. 
Deterministic  equations  have  been  developed  that  relate  a  media’s  permeability,  ep  and  attenuation2.  While  these 
equations  may  work  well  for  known  homogenous  media,  they  are  impossible  to  compute  for  dynamic  data  collections 
where  the  media’s  properties  are  unknown.  Instead  of  calculating  the  expected  attenuation  using  equations  dependent  on 
the  media’s  characteristics,  a  simple  attenuation  model  is  generated  directly  from  the  data.  This  model  is  obtained  by 
recording  the  maximum  amplitude  at  each  sample  time  across  all  traces  in  the  data  collection.  It  is  unlikely  that  this 
attenuation  model  will  have  the  monotonically  decreasing  curve  shape  required.  Smoothing  is  done  using  irregular  grid 
interpolation  between  the  maximum  monotonically  decreasing  peaks  of  the  curve  to  obtain  the  desired  shape.  The 
reciprocal  of  the  smoothed  data  is  used  to  create  a  time- varying  gain  function.  The  matched  filter  correlation  value  of  each 
detect  is  multiplied  by  the  value  of  this  gain  function  at  the  corresponding  sample  time  index  to  obtain  a  confidence  value. 
Finally,  only  the  detects  with  the  highest  confidence  values  are  kept,  depending  on  the  number  of  interfaces  desired. 

4.3  Prescreening 

To  characterize  layer  properties,  GPR  detections  must  be  segmented  according  to  the  layer  interfaces  that  caused  them. 
However,  to  improve  the  computational  efficiency  and  performance  of  the  segmentation/tracking  algorithms,  a  method  for 
rejecting  detects  that  do  not  belong  to  any  substantial  layer  interface  is  needed.  These  rejected  detects  are  considered  to  be 
clutter.  To  reject  clutter,  a  nonsupervised  classifier  is  used  to  determine  which  detects  do  not  have  properties  of  legitimate 
layers.  Classifier  performance  is  dependent  on  the  separability  of  classes  in  some  feature  space.  It  is,  therefore,  desirable 
to  obtain  a  set  of  features  for  each  candidate  detect  that  provides  some  separation  between  layers  and  clutter.  For 
simplicity  and  robustness,  two  easily  computable  features  are  used.  The  first  feature  considers  the  spatial  characteristics  of 
the  surrounding  detect  pattern  and  indicates  the  degree  to  which  the  shape  fits  that  of  a  layer.  The  second  feature  is  the 
detection  confidence  described  earlier.  The  extraction  of  the  layer  indication  features  and  the  classifier  procedures  are 
explained  in  this  section. 
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The  layer  indication  feature  is  a  combination  of  two  complementary  features.  Both  attempt  to  provide  some  general  cue  as 
to  the  probability  that  the  detect  of  interest  and  surrounding  detects  have  the  spatial  characteristics  of  a  layer.  The 
assumption  made  is  that  a  legitimate  layer  will  return  a  significant  density  of  detects  in  the  region  of  interest  which  have 
some  general  linear  dependence  with  regard  to  position.  The  two  features,  which  are  summed  to  create  the  layer 
indication  feature  for  classification,  are  a  histogram  indicator  and  a  correlation  indicator.  Uncorrelated  clutter  returns  have 
low  values  for  both  indicators 

The  histogram  indicator  feature  is  created  using  a  histogram  of  the  sample  index  of  all  the  detects  within  a  cross  range 
window  surrounding  the  detect  of  interest  as  illustrated  in  the  top  portion  of  Figure  3.  The  histogram  represents  the 
density  of  the  points  within  a  windowed  range  versus  sample  index.  The  density  value  at  the  sample  index  of  the  detect 
under  test  is  used  as  a  feature  called  the  histogram  indicator.  The  histogram  indicator  gives  a  high  value  for  artifacts  with 
well-defined  horizontal  characteristics,  such  as  layer  interfaces  with  minimal  slope. 

The  correlation  indicator  provides  better  indication  of  sloped-layer  interfaces.  It  considers  the  detects  within  a  window 
around  the  detect  under  test.  The  feature  value  used  is  the  absolute  value  of  the  correlation  between  the  down-range 
sample  index  and  the  cross-range  sample  index  of  the  group  of  windowed  detects  around  the  detect  under  test,  as  shown  in 
the  bottom  of  Figure  3.  Detects  from  sloped  layers  will  have  a  linear  dependence  between  these  indexes  and,  hence,  a  high 
correlation  value.  Returns  from  minimally-sloped  layers  tend  to  vary  randomly  around  the  mean  detect  index  of  the  layer, 
creating  a  slightly  lower  correlation  value  for  horizontal  layers.  Sloped  layers  typically  are  steep  enough  to  provide  a 
monotonic  relationship  between  their  indexes  in  spite  of  random  variance  in  detect  location  and  preserve  high  correlation 
values.  The  histogram  indicator  feature  and  the  correlation  indicator  feature  are  added  together  to  form  the  layer  indicator. 
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Figure  3.  Histogram  indicator  (top),  correlation  indicator  (bottom) 
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After  the  layer  indicator  feature  and  detection  confidence  are  obtained,  they  are  normalized  and  scaled  to  unity,  creating 
the  2-D  feature  space  for  separating  clutter  from  layer  detects.  The  features  are  supplied  to  a  K-means  clustering 
algorithm  which  autonomously  separates  the  clusters  based  on  the  minimization  of  a  performance  index8.  In  this  case,  the 
performance  index  is  defined  as  the  sum  of  the  squared  distances  from  all  points  in  a  cluster  domain  to  the  cluster  center 
K-means  clustering  segments  the  feature  space  to  allow  detects  that  have  clutter  characteristics  (low  values)  to  be 
discarded  prior  to  track  processing.  Essentially,  K-means  provides  an  autonomous  threshold  selection  based  on  the  natural 
distribution  of  the  detect  features.  The  nonsupervised  nature  of  this  algorithm  allows  it  to  operate  unmodified,  regardless 
of  the  GPR  system  or  collection  configuration  used. 

4.4  Tracking 

GPR  operators  can  often  associate  groups  of  detects  from  a  layer  interface  or  anomaly  by  visual  inspection.  However, 
computers  that  process  GPR  data  must  have  a  way  of  associating  layer  interfaces  and  geophysical  features  together  so  that 
physical  measurements  and  statistical  analysis  can  be  automated.  The  tracking  algorithm  assigns  a  unique  label  called  a 
track  number  to  groups  of  associated  detects.  The  tracking  algorithm  is  based  on  the  weighted  window  searching  and 
grouping  procedure  described  by  Cowdery9,  and  operates  on  the  detects  that  pass  the  prescreening  stage. 

After  tracking  is  completed,  there  are  often  several  separate  tracks  that  can  be  identified  as  belonging  to  one  consistent 
layer  interface.  Breaks  in  the  tracks  can  occur  due  to  misclassification  of  detects  by  the  clutter  rejection  algorithm  because 
tracking  assigns  track  labels  only  to  those  detects  that  are  not  classified  as  clutter.  Track  joining  works  by  considering  the 
clutter  detects  as  potential  layer  interfaces  and  assigning  track  labels  where  appropriate.  By  reassigning  labels  to  these 
misclassified  detects,  track  segments  can  be  automatically  joined  together.  The  joining  algorithm  starts  by  performing  an 
Euclidian  distance  search  for  nearby  detects  at  the  end  of  each  track  segment.  The  search  is  performed  inside  a  windowed 
area  determined  by  the  tracking  window  sizes.  If  a  non-tracked  detect  is  found,  it  is  assigned  to  the  current  track  and  the 
search  begins  again  at  the  newly  added  detect  location.  When  another  track  segment  is  found,  it  is  relabeled  with  the  track 
number  belonging  to  the  segment  where  the  search  was  initiated.  If  no  detects  or  tracks  are  found,  the  search  moves  to  the 
next  track  segment. 

Often,  returns  corresponding  to  legitimate  layers  are  discarded  by  the  detection  algorithm.  Track  filling  uses  a  windowed 
search  from  the  start  to  end  of  each  track  segment  and  attempts  to  fill  in  any  discontinuities  along  the  length  of  the  track 
segment  using  the  locations  of  peak  values  from  the  received  waveforms’  matched-filter  response.  Track  filling  and 
joining  increases  track  continuity  in  most  cases. 


5.  LAYER  PROPERTY  ANALYSIS 

Once  layers  have  been  segmented,  properties  of  the  layers  may  be  computed.  Of  primary  interest  is  the  dielectric  constant 
of  the  layer,  which  determines  signal  velocity  in  a  medium  and  is  combined  with  time-delay  measurements  to  an  interface 
for  calculation  of  interface  depth  and  layer  thickness.  Other  calculations  of  interest  may  include  statistical  characteristics 
of  layer  thickness  or  depth,  layer  length,  or  perhaps  even  surface  roughness  and  void  content  estimates.  Arriving  at  the 
dielectric  constant  of  layers  is  accomplished  in  different  ways  depending  upon  which  GPR  system  configuration  is 
utilized.  Two  methods/systems  are  used  here:  single-channel,  air-launched,  amplitude-based;  and  dual-channel,  ground- 
coupled,  time-based.  Each  system  presents  its  own  set  of  advantages,  disadvantages,  and  processing  requirements. 

5.1  Air-launched  system 

The  determination  of  layer  thickness  or  anomaly  location  is  a  multi-step  process  that  begins,  directly  or  indirectly,  with 
the  computation  of  the  reflection  coefficient,  p.  At  the  boundary  between  two  mediums  (layers,  voids,  etc.)  of  differing  ep 
p  is  the  ratio  of  the  reflected  to  incident  electric  fields.  Neglecting  multiple  reflections  at  each  boundary,  we  can  relate  p  at 
the  boundary  of  layer  i,  /+7  to  the  received  amplitudes  using  (2),  where  p0  =  -  A0  /  Am  and  Am  is  the  peak  amplitude  of  a 
flat  metal  plate  response.  The  dielectric  constant  of  layer  /+/  can  now  be  found  using  (3)  where  e0  =  dielectric  constant  of 
air  =  1 .  Finally,  the  layer  thickness  is  found  using  (4)  where  v,-  is  the  velocity  of  propagation  of  the  pulse  in  layer  i  and  c 
is  the  speed  of  light. 
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5.2  Ground-coupled  system 

Ground-coupled  systems  do  not  provide  an  accurate  surface  response  from  which  amplitude  measurements  may  be 
utilized  in  the  dielectric  calculations  for  successive  layers.  For  this  reason,  an  alternate  method  of  determining  the 
dielectric  constant  and  signal  velocity  through  each  layer  must  be  used.  Time  measurements  through  two  signal- 
propagation  paths  may  be  used  to  compute  the  average  velocity  to  a  layer  interface. 


H -  s  - H 


T 

d 
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Figure  4.  Dual  antenna  configuration 

For  a  fixed  geometry  configuration,  as  shown  in  Figure  4,  the  separation  distance  between  the  antennas,  s ,  and  the  travel 
time  through  each  path  are  known,  leaving  only  distance  and  velocity  to  each  interface  as  unknowns.  The  direct  path 
distance  r  and  cross  path  distance  d  are  given  by  (5)  and  (6),  where  tr  and  t#  are  the  measured  2-way  path  propagation 
times.  Both  distance  terms  can  be  solved  for  in  the  two-equation  system  with  the  average  velocity  v  calculated  using  (7). 
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,  where  v  is  the  average  velovity. 
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In  applying  this  technique,  subsurface  profiles  are  collected  through  multiple  propagation  paths  at  each  sample  location. 
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Each  channel  of  profile  data  is  then  independently  processed  by  the  detection  and  tracking  algorithms.  This  results  in 
detected  interfaces  in  each  profile.  The  time  to  an  interface  in  each  channel  will  be  used  to  determine  propagation  velocity 
and  depth;  however,  up  to  this  point,  no  associations  have  been  made  between  detects  from  an  interface  in  one  antenna 
channel  (or  path)  and  the  detects  from  the  same  interface  in  the  second  antenna  channel.  It  is,  therefore,  critical  that  some 
method  pair  together  corresponding  detects  in  profile  data  from  each  channel.  Wave  velocity,  medium  dielectric  constant, 
and  interface  depth  may  then  be  calculated. 

5.2.1  Matching  tracks 

Matching  is  a  processing  stage  that  is  unique  to  the  dual-antenna  system.  Layer  dielectric  constant  is  determined  by 
utilizing  time  measurements  through  two  wave  propagation  paths  in  the  two-antenna  system.  This  requires  detected  layers 
in  the  data  from  one  antenna  to  be  paired  together  with  the  corresponding  detected  layer  in  the  data  from  the  second 
antenna.  An  algorithm  has  been  designed  which  autonomously  performs  this  pairing  based  upon  acceptable  locations 
according  to  a  range  of  possible  dielectric  values. 

Calculations,  given  the  system  geometry,  were  performed  to  determine  the  range  of  acceptable  delays  between  channels 
versus  direct-channel  sample  time  for  an  assumed  range  of  dielectric  values.  From  this  range  of  acceptable  delays 
between  channels,  a  minimum  cross-channel  offset  curve  was  determined.  This  curve  provides  the  offset  position  of  the 
matching  window  in  the  cross  channel  with  respect  to  the  sample  position  in  the  direct  channel.  The  matching  window 
size  is  also  determined  from  the  range  of  acceptable  delays.  The  positions  of  tracked  detects  in  the  direct  versus  cross 
channels  are  compared  to  determine  matched  pairs.  The  algorithm  biases  matching  towards  the  lowest  er  selection. 

5.2.2  Depth  calculation 

To  calculate  the  depth  of  each  matched  pair  of  detects,  the  detect  location  index  in  the  two  channels  are  used  to  reference 
a  calibrated  matrix  of  average  e,.  values,  which  are  unique  for  a  given  pair  of  channel  propagation  times.  The  calibrated 
matrix  was  generated  from  testbed  measurement  calculations  using  materials  with  known  e,.  over  a  range  of  depth 
intervals.  The  average  ^obtained  from  the  calibrated  matrix  is  then  used  to  determine  the  average  velocity  of  the  radar 
signal  to  the  interface  of  interest.  From  this  velocity  and  the  time  through  either  channel,  the  depth  to  the  interface  is 
computed. 

5.2.3  Layer  dielectric  constant  calculation 

A  parameter  of  interest  is  the  dielectric  constant,  Ep  of  each  detected  layer.  Only  the  average  dielectric  constant  to  each 
interface  has  been  computed  (to  this  point  in  the  signal  processing  chain)  because  the  time  to  each  interface,  as  measured 
by  the  GPR,  is  a  function  of  the  average  velocity  of  wave  propagation  to  the  interface  of  interest.  The  dielectric  of  the 
medium  above  an  interface  is,  however,  often  more  useful  information  than  the  average  dielectric  of  all  mediums  above 
the  interface.  For  this  reason,  an  algorithm  has  been  written  that  uses  an  iterative  method,  starting  with  the  first  layer  and 
working  down,  to  calculate  an  estimate  of  the  dielectric  of  each  individual  layer.  The  algorithm  assumes  that  the  first 
matched  pair  in  a  given  trace  (a  single  GPR  measurement  of  the  subsurface  at  one  Iocation-sometimes  called  a  profile; 
although,  a  profile  may  sometimes  refer  to  a  succession  of  traces  stacked  together  to  represent  an  area  of  the  subsurface) 
is  from  the  first-to-second  layer  interface;  thus,  the  average  dielectric  constant  is  equal  to  the  actual  dielectric  of  the  first 
layer.  Once  the  first  layer  parameters  have  been  established,  the  next  layer  (matched  pair)  is  addressed.  Given  the  depth 
and  average  dielectric  to  each  interface  and  preceding  layer,  the  dielectric  constant  of  the  layer  may  be  calculated  using 
(8)  and  (9)  where  di ,  W,  and  are  the  depth,  velocity,  and  er  respectively  for  each  layer  i.  This  process  is  successively 
followed  for  each  interface  in  a  trace. 
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6.  DATA  ANALYSIS  EXAMPLES 


A  pair  of  GPR  data  collections  using  the  ground-coupled  configuration  described  earlier  were  recently  completed  in 
conjunction  with  two  traditional  airfield  pavement  evaluation  surveys  performed  by  AFSECA  personnel.  These 
collections  provided  a  large  set  of  ground-truthed  GPR  data  needed  for  algorithm  development,  verification,  and  fine 
tuning.  The  large  number  of  data  files  and  other  associated  survey  data  obtained  also  prompted  software  interface 
adaptations  to  the  GPR  analysis  software  to  provide  more  efficient  processing,  analysis  and  interaction  with  the  large 
number  of  data  sets.  The  figures  below  show  data  from  one  of  these  collections,  along  with  the  detected  interface 
boundaries  and  plots  of  the  calculated  depths.  The  GPR  analysis  software  was  used  to  process  and  display  the  data, 
remove  undesired  detected  interfaces,  and  generate  the  figures  shown  below  in  Figure  5. 


Figure  5.  GPR  data  from  airfield  at  McChord  AFB.  The  two  figures  on  the  left  show  GPR  data  collected  from  a  pair  of  ground-coupled 
antennas.  The  image  on  the  top  is  from  the  direct  channel  and  the  image  on  the  bottom  is  from  the  cross  channel.  The  two  figures  in  the 
middle  show  the  same  data  sets  with  the  first  two  detected  layer  interface  boundaries  over-plotted.  The  plots  on  the  right  show  the 
calculated  depth  of  the  first  two  interface  boundaries  (note  that  due  to  the  method  used  to  calculate  depth,  the  plots  are  identical).  The 
black  line  and  arrow  in  each  of  the  figures  represents  the  closest  point  in  the  data  to  a  collected  core  sample.  For  this  example,  the 
thickness  of  the  pavement  core  was  5.5  inches. 

Examples  of  the  GPR  pavement  thickness  calculations  versus  nearby  measured  cores  from  the  airfield  evaluation  at 
McChord  AFB  are  shown  in  Table  1.  The  results  were  obtained  after  processing  the  GPR  data  using  the  algorithms 
described  in  this  paper.  The  GPR  analysis  software  was  used  to  remove  undesired  detected  interfaces  and  connecting  non- 
continuous  tracked  segments,  as  necessary.  For  some  data  sets,  layer  interfaces  can  remain  undetected  due  to  either  weak 
signal  response  (core  3)  or  interference  from  multiple  thin  layers  (core  10).  When  the  pavement  interface  can  be  detected, 
the  calculated  thickness  values  were  close  to  core  measurements. 
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Table  1.  Measured  core  thickness  vs.  calculated  GPR  thickness  from  airfield  at  McChord  AFB 


Core 

Number 

Measured 
Thickness  (In.) 
(nearest  .25  In.) 

GPR  Calculated 
Thickness  (In.) 

Pavement 

Type 

Notes 

i 

12.0 

11.87 

RPCC 

Rebar  noted  3.5  In.  from  top 

•* 

2 

12.0 

11.07  ■ 

RPCC 

Rebar  noted  3.5  In.  from  to 

3 

11.50 

N/A 

PCC 

Pavement  layer  not  detected  in  cross  channel. 

4 

11.75 

11.74 

PCC 

5 

5.75 

6.65 

AC 

6 

5.50 

5.32 

AC 

7 

5.50 

5.20 

AC 

8 

5.25 

5.29 

AC 

9 

5.25 

5.38 

AC 

10 

5.25 

N/A 

AC 

Pavement  layer  not  detected  in  cross  channel. 

RPCC  -  Reinforced  Portland  Cement  Concrete,  PCC  -  Portland  Cement  Concrete,  AC  -  Asphaltic  Concrete 

7.  SUMMARY 

GPR  has  the  potential  the  improve  airfield  evaluations  by  allowing  more  continuous  and  rapid  assessment  of  pavement 
and  subsurface  properties.  The  signal  processing  methods  described  in  this  paper  provide  estimates  of  the  properties  of 
pavement  and  other  subsurface  layers.  Large  amounts  of  GPR  data  can  be  quickly  processed  and  the  results  can  be  easily 
analyzed  and  summarized  using  the  GPR  analysis  software.  When  favorable  GPR  conditions  exist,  the  estimated 
pavement  thickness  is  generally  in  close  agreement  with  the  measured  values. 
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